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We study the dynamics of confinement-deconfinement (C-D) phase transition in the context of 

relativistic heavy-ion colhsions within the framework of effective models for the Polyakov loop order 

parameter. We study the formation of Z{3) walls and associated strings in the initial transition from 

the confining (hadronic) phase to the deconfining (QGP) phase via the so called Kibble mechanism. 

Essential physics of the Kibble mechanism is contained in a sort of domain structure arising after 

f^ any phase transition which represents random variation of the order parameter at distances beyond 

T-H the typical correlation length. We implement this domain structure by using the Polyakov loop 

^^ effective model with a first order phase transition and confine ourselves with temperature/time 

Cn ranges so that the first order C-D transition proceeds via bubble nucleation, leading to a well 

,—1 defined domain structure. The formation of Z{3) walls and associated strings results from the 

3 coalescence of QGP bubbles expanding in the confining background. We investigate the evolution 

' of the Z{3) wall and string network. We also calculate the energy density fiuctuations associated 

00 with Z{3) wall network and strings which decay away after the temperature drops below the quark- 

CNj hadron transition temperature during the expansion of QGP. We discuss evolution of these quantities 

with changing temperature via Bjorken's hydrodynamical model and discuss possible experimental 

' —I ' signatures resulting from the presence of Z{3) wall network and associate strings. 
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^D Search for the quark-gluon plasma (QGP) at relativistic heavy-ion collision experiments (RHICE) has reached a 

^^ very exciting stage with the ongoing experiments at RHIC and with the upcoming heavy-ion experiments at the Large 
Hadron Collider. No one doubts that the QGP phase has already been created at RHIC but conclusive evidence for 
t^ the same is still lacking. Many signatures have been proposed for the detection of the QGP phase [T] and these 
^P have been thoroughly investigated both theoretically and experimentally. Along with continued investigation of these 
important signatures of QGP, there is a need for investigating novel signals exploring qualitatively non-trivial features 
of the QGP phase and/or the quark-hadron phase transition. 

With this view we focus on the non-trivial vacuum structure of the QGP phase which arises when one uses the 
expectation value of the Polyakov loop l{x) as the order parameter for the confinement-deconfinement phase transition 
[2]. This order parameter transforms non-trivially under the center Z{3) of the color SU(3) group and is non-zero 
C^ above the critical temperature Tc. This breaks the global Z(3) symmetry spontaneously above Tc, while the symmetry 

is restored below Tc in the confining phase where this order parameter vanishes. In the QGP phase, due to spontaneous 
breaking of the discrete Z{3) symmetry, one gets domain walls (interfaces) which interpolate between different Z{3) 
vacua. The properties and physical consequences of these Z(3) interfaces have been discussed in the literature 3J. It 
has been suggested that these interfaces should not be taken as physical objects in the Minkowski space [2|. Similarly, 
it has also been subject of discussion whether it makes sense to talk about this Z(3) symmetry in the presence of 
quarks [5 . However, we will follow the approach where the presence of quarks is interpreted as leading to explicit 
breaking of Z{3) symmetry, lifting the degeneracy of different Z{3) vacua [6}jE]. Thus, with quarks, even planar Z(3) 
interfaces do not remain static and move away from the region with the unique true vacuum. Our main discussion 
will be for the pure gauge theory which we discuss first. Later we will briefly comment on the situation with quarks. 
A detailed study with inclusion of quark effects is postponed for a future work. 
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In earlier works some of us had shown that there are novel topological string defects which form at the intersection 
of the three Z{3) interfaces. These strings are embedded in the QGP phase and their cores consist of the confining 
phase. Structure of these strings and interfaces were discussed in these earlier works (TUI. It was also shown that 
reflection of quarks from collapsing Z{3) interfaces can lead to large scale baryon inhomogeneities in the early universe 
[llj . This effect is also utilized to argue for Pt enhancement of quarks of heavy flavor (consecutively corresponding 
hadrons) in relativistic heavy-ion collisions [12) . 

In this paper we carry out numerical simulation of formation of these Z{3) interfaces and associated strings at the 
initial confinement-deconfinement transition which is believed to occur during the pre-equilibrium stage in relativistic 
heavy- ion collision experiments. For the purpose of numerical simulation we will model this stage as a quasi-equilibrium 
stage with an effective temperature which first rises (with rapid particle production) to a maximum temperature 
Tq > Tc, where T^. is the critical temperature for the confinement-deconfinement phase transition, and then decreases 
due to continued plasma expansion. 

We use the effective potential for the Polyakov loop expectation value l{x) as proposed by Pisarski [71|H] to study 
the confinement-deconfinement (C-D) phase transition. Within this model, the C-D transition is weakly first order. 
Even though Lattice results show that quark-hadron transition is most likely a smooth cross-over at zero chemical 
potential, in the present work we will use this first order transition model to discuss the dynamical details of quark- 
hadron transition. One reason for this is that our study is in the context of relativistic heavy-ion collision experiments 
(RHICE) where the baryon chemical potential is not zero. For not too small values of the chemical potential, the 
quark-hadron phase transition is expected to be of first order, so this may be the case relevant for us any way 
(especially when collision energy is not too high). Further, our main interest is in determining the structure of the 
network of Z(3) domain walls and strings resulting during the phase transition. These objects will form irrespective of 
the nature of the transition, resulting entirely from the finite correlation lengths in a fast evolving system, as shown 
by Kibble [13] . The Kibble mechanism was first proposed for the formation of topological defects in the context of 
the early universe [13J, but is now utilized extensively for discussing topological defects production in a wide variety 
of systems from condensed matter physics to cosmology [Mj. Essential ingredient of the Kibble mechanism is the 
existence of uncorrelated domains of the order parameter which result after every phase transition occurring in finite 
time due to finite correlation length. A first order transition allows easy implementation of the resulting domain 
structure especially when the transition proceeds via bubble nucleation. With this view, we use the Polyakov loop 
model as in T, "8; to model the phase transition and confine ourselves with temperature/time ranges so that the first 
order quark-hadron transition proceeds via bubble nucleation. 

The Z{3) wall network and associated strings (as mentioned above) formed during this early confinement- 
deconfinement phase transition evolve in an expanding plasma with decreasing temperature. Eventually when the 
temperature drops below the deconfinement-confinement phase transition temperature Tc, these Z{3) walls and associ- 
ated strings will melt away. However they may leave their signatures in the form of extended regions of energy density 
fiuctuations (as well as Pt enhancement of heavy-flavor hadrons [T2]). We make estimates of these energy density 
fluctuations which can be compared with the experimental data. Especially interesting will be to look for extended 
regions of large energy densities in space-time reconstruction of hadron density (using hydrodynamic models) . In our 
model, we expect energy density fluctuations in event averages (representing high energy density regions of domain 
walls/strings), as well as event-by-event fiuctuations as the number/geometry of domain walls/strings and even the 
number of QGP bubbles, varies from one event to the other. A detailed analysis of energy fiuctuations, especially the 
event-by-evcnt fiuctuations, is postponed for a future work. 

We also determine the distribution/shape of Z{3) wall network and its evolution. In particular, our results provide 
an estimate of domain wall velocities (for the situations studied) to range from 0.5 to 0.8. These results provide crucial 
ingredients for a detailed study of the effects of collapsing Z{3) walls on the Pt enhancement of heavy flavor hadrons 
[TTIIT2] in RHICE. We emphasize that the presence of Z{3) walls and string may not only provide a qualitatively new 
signature for the QGP phase in these experiments, it may provide the flrst (and may be the only possible) laboratory 
study of such topological objects in a relativistic quantum fleld theory system. 

The paper is organized in the following manner. In section II, we discuss the Polyakov loop model of confinement- 
deconfinement phase transition. We describe the effective potential proposed by Pisarski 7J and discuss the structure 
of the Z{3) walls and associated strings. In section HI, we present the physical picture of the formation of these 
Z(3) walls and associated strings in the confinement-deconfinement phase transition via the Kibble mechanism which 
provides a general framework for the production of topological defects in symmetry breaking phase transitions. We 
confine ourselves with temperature/time ranges so that the first order transition (in the Pisarski model) proceeds 
via bubble nucleation. The other possibility of spinodal decomposition is of completely different nature and we will 
present it in a future work. (We mention here that a simulation of spinodal decomposition in Polyakov loop model 
has been carried out in ref. [15] . where fluctuations in the Polyakov loop are investigated in detail. In comparison, the 
main focus of our work is on the topological objects, Z(3) walls, and strings, and energy density fluctuations resulting 
therefrom.) In section IV, we discuss the calculation of the proflle of the critical bubble using bounce technique [T5] 



and also the estimates of nucleation rates for bubbles for the temperature/time range relevant for RHICE. Section V 
presents the numerical technique of simulating the phase transition via random nucleation of bubbles. In section VI 
we discuss the issue of the effects of quarks in our model. Section VII presents the results of the numerical simulations. 
Here we discuss distribution of Z{3) wall/string network formed due to coalescence of QGP bubbles (in a confining 
background). We also calculate the energy density fluctuations associated with Z{3) wall network and strings. We 
discuss evolution of these quantities with changing temperature via Bjorken's hydrodynamical model. In section VIII 
we discuss possible experimental signatures resulting from the presence of Z(S) wall network and associate strings. 
Section IX presents conclusions. 

II. THE POLYAKOV LOOP MODEL 

We first discuss the case of pure SU(N) gauge theory. We will later discuss briefly the case with quarks. The order 
parameter for the confinement-deconfincincnt (C-D) phase transition is the expectation value of the Polyakov loop 
l{x) which is defined as 



l{x) = {l/N)tr{Pexp{ig Ao{x,T)dT)) (1) 

Jo 

Here P denotes the path ordering, g is the gauge coupling, and (3 = l/T, where T is the temperature. Aq{x,t) = 
^g(x,r)ra, Ta being the generators of SU(N) in the fundamental representation, is the time component of the vector 
potential at spatial position x and Euclidean time r. Under a global Z(N) symmetry transformation, l{x) transforms 
as 

l{x) -> exp{2TTm/N)l{x), n :^ 0, 1 , ...{N - 1) (2) 

The expectation value of l{x) is related to e~^^ where F is the free energy of an infinitely heavy test quark. For 
temperatures below Tc, in the confined phase, the expectation value of Polyakov loop is zero corresponding to the 
infinite free energy of an isolated test quark. (Hereafter, we will use the same notation l{x) to denote the expectation 
value of the Polyakov loop.) Hence the Z(N) symmetry is restored below Tc. Z(N) symmetry is broken spontaneously 
above Tc where l{x) is non-zero corresponding to the finite free energy of the test quark. For QCD, N = 3, and we 
take the effective theory for the Polyakov loop as proposed by Pisarski [71 18] . The effective Lagrangian density is 
given by 

L^^\d,l\'T^-Vil) (3) 

V{1) is the effective potential for the Polyakov loop 

v{i) = {^\i\' - ^f{i' + nr) + ^{\i\'r)b4 T^ (4) 

The values of various parameters are fixed to reproduce the lattice results [T71 [TB] for pressure and energy density 
of pure SU(3) gauge theory. We make the same choice and give those values below. The coefficients 63 and 64 have 
been taken as, 63 — 2.0 and 64 — 0.6016. We will take the same value of 62 for real QCD (with three massless quark 
flavors), while the value of 64 will be rescaled by a factor of 47.5/16 to account for the extra degrees of freedom relative 
to the degrees of freedom of pure SU(3) gauge theory [T7]. The coefficient 62 is temperature dependent pTj. 62 is 
taken as, 62(r) = (1- l.ll/r)(l + 0.265/r)2(l + 0.3/r)3-0.487, where r is taken as T/T^. With the coefficients chosen 
as above, the expectation value of order parameter approaches to x = ^3/^ + 4 v ^| + ^&2(T — 00) for temperature 
T — > 00. As in fS], we use the normalization such that the expectation value of order parameter Ig goes to unity 
for temperature T — ?> 00. Hence the fields and the coefficients in V{1) are rescaled as I — ^ l/x, 62(2^) — ^ b2{T)/x'^, 
^3 -^ b^/x and 64 — > b^x'^ to get proper normalization of Ig. 

For the parameters chosen as above, the value of Tc is taken to be 182 MeV. We see that the 63 term in Eq.(4) gives 
a cos(36') term, leading to Z(3) degenerate vacua structure. Here, the shape of the potential is such that there exists 
a metastable vacuum upto a temperature ^ 250 MeV. Hence first order transition via bubble nucleation is possible 
only upto T = 250 MeV. We show the plot of V{1) in 6* = direction in Fig. (la) for a value of temperature T = 
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FIG. 1: 

(a) Plot of V{V) in S = direction for T = 185AfeV showing the metastable vacuum at / = 0. In (a) and (b), plots of V are 

given in units of T^ . The value of critical temperature is taken to be Tc — 182MeV. The Z(S) structure of the vacuum can be 

seen in (b) in the plot of the potential V{1) as a function of 9 for fixed j/| — Iq. Here, lo corresponds to the absolute minimum 

of V{1). 



185 MeV. This shows the metastable vacuum at I = 0. Fig. (lb) shows the structure of vacuum by plotting V(/) as a 
function of 6 for fixed Ig, where ^o is the vacuum expectation value of V{1) at T — 185 MeV. 

In Fig. (lb), the three degenerate vacua are separated by large barrier in between them. While going from one 
vacuum to another vacuum, the field configuration is determined from the field equations. For the case of degenerate 
vacua, there are time independent solutions which have planar symmetry. These solutions are called domain wall. For 
the non-degenerate case, as will be appropriate for the case when quarks are included as dynamical degrees of freedom 
in discussing the quark hadron transition, the solutions of the interfaces separating these vacua will be similar to 
the bounce solutions \W , though the standard bounce techniques need to be extended for the case of complex scalar 
field. The resulting planar domain wall solutions will not be static. As mentioned above, we will be neglecting such 
effects of quarks, and hence will discuss the case of degenerate vacua only. Later, in section VI, we will examine the 
justification of using the approximation of neglecting quark effects. 

In physical space, after the phase transition, regions with different Z(3) vacua are separated by domain walls. Inside 
a domain wall, \l\ becomes very small as T — >■ T^ [10]. In ref.[TD] the intersection of the three Z(3) domain walls was 
considered and, using topological arguments, it was then shown that at the line-like intersection of these interfaces, 
the order parameter l{x) should vanish. This leads to a topological string configuration with core of the string being 



in the confining phase. Properties of these new string configurations were determined in [TO] using the model of 
Eq.(3). This string configuration is interesting, especially when we note that such a string has exactly reverse physical 
behavior compared to the standard QCD string. The QCD string exists in the confining phase, connecting quarks 
and antiquarks, or forming baryons, glueballs etc. Inside the QCD string, the core region is expected to behave as a 
deconfined region. In contrast the string discussed in 10:, arising at the intersection of Z(3) walls, exists in the high 
temperature deconfined phase. Its core is characterized by restored Z{3) symmetry, implying that it is in the confined 
phase. To differentiate it with the standard QCD string, this new string structure was called as the QGP string in 
ref . |10) . It is also important to note that although the standard QCD string breaks by creating quark-antiquark pairs, 
the QGP string cannot break as it originates from topological arguments. This QGP string, thus should either form 
closed loops, or it should end at the boundary separating the deconfined phase from the confined phase. Its structure is 
very similar to certain axionic strings discussed in the context of early universe |19| . Note that as these strings contain 
the confining phase (with I = 0) in the core, while they are embedded in the QGP phase, a transition from deconfining 
phase to the confining phase, in the presence of such strings may begin from regions near the strings. Similarly, the 
presence of domain walls may lead to heterogeneous bubble nucleation in a first order quark-hadron transition. We 
will be studying the formation of Z(3) domain walls and these QGP strings in the initial confinement-deconfinement 
phase transition, and their subsequent evolution. 

III. DOMAIN WALL AND STRING FORMATION VIA KIBBLE MECHANISM 

We now briefly describe the physical picture of the formation of topological defects via the Kibble mechanism. 
Kibble first gave a detailed theory of the formation of topological defects in symmetry breaking phase transitions 
in the context of the early universe [TOj. Subsequently it was realized that the basic physics of Kibble mechanism 
is applicable to every symmetry breaking transition, from low energy physics of condensed matter systems to high 
energy physics relevant for the early universe [Tl ■ 

Basic physics of the Kibble mechanism can be described as follows. After a spontaneous symmetry breaking phase 
transition, the physical space consists of regions, called domains. In each domain the configuration of the order 
parameter field can be taken as nearly uniform while it varies randomly from one domain to another. In a numerical 
simulation where the phase transition is modeled to implement the Kibble mechanism, typically the physical region 
is divided in terms of elementary domains of definite geometrical shape. The order parameter is taken to be uniform 
within the domain and random variations of order parameter field within the vacuum manifold are allowed from one 
domain to the other. The order parameter field configuration in between domains is assumed to be such that the 
variation of the order parameter field is minimum on the vacuum manifold (the so-called geodesic rule) . With this 
simple construction, topological defects arise at the junctions of several domains if the variation of the order parameter 
in those domains traces a topologically non-trivial configuration in the vacuum manifold. See, ref. j20| for a detailed 
discussion of this approach. 

We, however, will follow a more detailed simulation as in ref.j21| where the Kibble mechanism was implemented in 
the context of a first order transition. Bubbles of true vacuum (determined from the bounce solution) were randomly 
nucleated in the background of false vacuum. Each bubble was taken to have the uniform orientation of the order 
parameter in the vacuum manifold, while the order parameter orientation varied randomly from one bubble to another. 
This provided the initial seed domains, as needed for the Kibble mechanism. Evolution of this initial configuration 
via the field equation led to expansion of bubbles which eventually coalesce and lead to the formation of topological 
defects at the junctions of bubbles when the order parameter develops appropriate variation (winding) in that region. 
Important thing is that in this case one does not need to assume anything like the geodesic rule. As different bubbles 
come into contact during their expansion, the value of the order parameter in the intermediate region is automatically 
determined by the field equations. 

An important aspect of the Kibble mechanism is that it does not crucially depend on the dynamical details of the 
phase transition. Although the domain size depends on the dynamics of phase transitions, the defect number density 
(per domain) and type of topological defects produced via the Kibble mechanism depends only on the topology of order 
parameter space and spatial dimensions. If the vacuum manifold M has disconnected components, then domain walls 
form. If it is multiply connected (i.e., if M contains unshrinkable loops), strings will form. When M contains closed 
two surfaces which cannot be shrunk to a point, then monopoles will form in three dimensional physical space. In our 
case domain walls and string network will be produced in the QGP phase. As discussed above, domain walls arise 
due to interpolation of field between different Z(3) vacua. At the intersection of these interfaces, string is produced. 



IV. CRITICAL BUBBLE PROFILE AND NUCLEATION PROBABILITY 

As we explained in the introduction, we will study the Z(3) domain wall and string formation with the first order 
transition model given by Eq.(3), such that the transition occurs via bubble nucleation. The semiclassical theory of 
decay of false vacuum at zero temperature has been given in ref. |16| . and the finite temperature extension of this 
theory was given in ref. |22) . The process of barrier tunneling leads to the appearance of bubbles of the new phase. 
Resulting bubble profile is determined using the bounce solution for the false vacuum decay which we will discuss 
below. 

First we note general features of the dynamics of a standard first order phase transition at finite temperature via 
bubble nucleation. A region of true vacuum, in the form of a spherical bubble, appears in the background of false 
vacuum. The creation of bubble leads to the change in the free energy of the system as, 

A'jr 

F{R)^ F^ + F^^A-KR^a- —R^T] (5) 

o 

where R is the radius of bubble, Fs is the surface energy contribution and F„ is the volume energy contribution, rj 
is the difference of the free energy between the false vacuum and the true vacuum and a is the surface tension which 
can be determined from the bounce solution. 

A bubble of size R will expand or shrink depending on which process leads to lowering of the free energy given 
above. The bubbles of very small sizes will shrink to nothing since surface energy dominates. If the radius of bubble 
exceeds the critical size Re ~ — , it will expand and lead to the transformation of the mctastable phase into the stable 
phase. 

Eq.(5) is useful for the so-called thin wa// bubbles where there is a clear distinction between the surface contribution 
to the free energy and the volume contribution. For the temperature/time relevant for our case in relativistic heavy-ion 
collisions this will not be the case. Instead we will be dealing with the thick wall bubbles where surface and volume 
contributions do not have clear separations. We will determine the profiles of these thick wall bubbles numerically 
following the bounce technique [16J. 

First we note that for the effective potential in Eq. (4) , the barrier between true vacuum and false vacuum vanishes 
at temperatures above about 250 MeV. So first order transition via bubble nucleation is possible only within the 
temperature range of T^. ~ 182 MeV - 250 MeV. Above T ~ 250 MeV, spinodal decomposition will take place due to 
the roll down of field. Implementation of the dynamics of phase transition via spinodal decomposition is of completely 
different nature, and we hope to discuss this in a future work. 

As we are discussing the initial confinement-deconfinement transition in the context of RHICE, clearly the discussion 
has to be within the context of longitudinal expansion only, with negligible effects of the transverse expansion. 
However, Bjorken's longitudinal scaling model [23 cannot be applied during this pre-equilibrium phase, even with the 
assumption of quasi-equilibrium (as discussed above), unless one includes a heat source which could account for the 
increase of effective temperature during this phase to the maximum equilibrium temperature Tq. As indicated above, 
this heat source can be thought of as representing the rapid particle production (with subsequent thermalization) 
during this early phase. We will not attempt to model such a source here. Instead, we will simply use the field 
equations resulting from Bjorken's longitudinal scaling model for the evolution of the field configuration for the entire 
simulation, including the initial pre-equilibrium phase from r = to t = tq. The heating of the system until t — tq 
will be represented by the increase of the temperature upto T = Tq. Thus, during this period, the energy density and 
temperature evolution will not obey the Bjorken scaling equations [53]. After tq, with complete equilibrium of the 
system, the temperature will decrease according to the equations in the Bjorken's longitudinal scaling model. 

We will take the longitudinal expansion to only represent the fact that whatever bubbles will be nucleated, they get 
stretched into ellipsoidal, and eventually cylindrical, shapes during the longitudinal expansion (ignoring the boundary 
effects in the longitudinal direction). The transverse expansion of the bubble should then proceed according to relative 
pressure difference between the false vacuum and the true vacuum as in the usual theory of first order phase transition. 
We will neglect transverse expansion for the system, and focus on the mid rapidity region. With this picture in mind, 
we will work with effective 2-f 1 dimensional evolution of the field configuration, (neglecting the transverse expansion 
of QGP). However, for determining the bubble profile and the nucleation probability of bubbles, one must consider 
full 3-1-1 dimensional case as bubbles are nucleated with full 3-dimensional profiles in the physical space. It will turn 
out that the bubbles will have sizes of about 1 - 1.5 fm radius. Taking the initial collision region during the pre- 
equilibrium phase also to be of the order of 1-2 fm in the longitudinal direction, it looks plausible that the nucleation 
of 3-dimensional bubble profile as discussed above may provide a good approximation. Of course the correct thing 
will be to consider the bounce solutions for rapidly longitudinally expanding plasma, and we hope to return to this 
issue in some future work. 



We neglect transverse expansion in the present work, which is a good approximation for the early stages when 
wall/string network forms. However, this will not be a valid approximation for later stages, especially when temper- 
ature drops below T^ and wall/string network melts. The way to account for the transverse expansion in the context 
of our simulation will be to take a lattice with much larger physical size than the initial QGP system size, and allow 
free boundary conditions for the field evolution at the QGP system boundary (which will still be deep inside the 
whole lattice). This will allow the freedom for the system to expand in the transverse direction automatically. With 
a suitable prescription of determining temperature from local energy density (with appropriate account of field con- 
tributions and expected contribution from a plasma of quarks and gluons) in a self consistent manner, the transverse 
expansion can be accounted for in this simulation. We hope to come back to this in a future work. 

Let us consider the effective potential in Eq.(4), at a temperature such that there is a barrier between the true 
vacuum and the three Z(3) vacua. An example of this situation is shown in Fig.l for the case with T = 185 
MeV. The initial system (of nucleons) was at zero temperature with the order parameter l(x) = 0, and will be 
superheated as the temperature rises above the critical temperature. It can then tunnel through the barrier to 
the true vacuum, representing the deconfined QGP phase. At zero temperature, the tunneling probability can be 
calculated by finding the bounce solution which is a solution of the 4-dimensional Euclidean equations of motion 
However, at finite temperature, this 4-dimensional theory will reduce to an effectively 3 Euclidean dimensional theory 
if the temperature is sufficiently high, which we will take to be the case. 

For this finite temperature case, the tunneling probability per unit volume per unit time in the high temperature 
approximation is given by j22j (in natural units) 

T^A e-^3(/)/T (g) 

where 83(1) is the 3-dimensional Euclidean action for the Polyakov loop field configuration that satisfies the classical 
Euclidean equations of motion. The condition for the high temperature approximation to be valid is that T » r^ , 
where Tq is the radius of the critical bubble in 3 dimensional Euclidean space. The values of temperature for our 
case (relevant for bubble nucleation) will be above T = Tc = 182 MeV. As we will see, the bubble radius will be 
larger than 1.5 fm ( ~ (130 MeV)~^) which justifies our use of high temperature approximation to some extent. 
The determination of the pre-exponential factor is a non-trivial issue and we will discuss it below. The dominant 
contribution to the exponential term in F comes from the least action 0(3) symmetric configuration which is a solution 
of the following equation (for the Lagrangian in Eq.(3)). 

^-H^^^^^nO (7) 

dr^ rdr 2NT^ ^ ' ^ ' 



where r = te = \/ S?- + t\ , subscript E denoting the coordinates in the Euclidean space. 
The boundary conditions imposed on I are 

1 = as r — > 00 and -— = at r = (8) 

dr 

Bounce solution of Eq.(7) can be analytically obtained in the thin wall limit where the difference in the false vacuum 
and the true vacuum energy is much smaller than the barrier height. This situation will occur for very short time 
duration near T = Tc for the effective potential in Eq.(4). However, as the temperature is rapidly evolving in the 
case of RHICE, there will not be enough time for nucleating such large bubbles (which also have very low nucleation 
rates due to having large action). Thus the case relevant for us is that of thick wall bubbles whose profile has to be 
obtained by numerically solving Eq.(7). 

As we have mentioned earlier, in the high temperature approximation the theory effectively becomes 3 (Euclidean) 
dimensional. For a theory with one real scalar field in three Euclidean dimensions the pre-exponential factor arising 
in the nucleation rate of critical bubbles has been estimated, see ref. [22 j . The pre-exponential factor obtained from 
[22] for our case becomes 

It is important to note here that the results of [52] were for a single real scalar field and one of the crucial ingredients 
used in [32] for calculating the pre-exponential factor was the fact that for a bounce solution the only light modes 



contributing to the determinant of fluctuations were the deformations of the bubble perimeter. Even though we are 
discussing the case of a complex scalar field l{x), this assumption may still hold as we are calculating the tunneling 
from the false vacuum to one of the Z(3) vacua (which are taken to be degenerate here as we discussed above). This 
assumption may need to be revised when light modes e.g. Goldstone bosons are present which then also have to be 
accounted for in the calculation of the determinant. 

A somewhat different approach for the pre-exponential factor in Eq.(6) is obtained from the nuclcation rate of 
bubbles per unit volume for a liquid-gas phase transition as given in ref . |241 I25j , 

r = ^n,e-^^/^ (10) 

Here k is the dynamical prefactor which determines the exponential growth rate of critical droplets, fig is a 
statistical prefactor which measures the available phase space volume. The exponential term is the same as in Eq. (6) 
with AF being the change in the free energy of the system due to the formation of critical droplet. This is the 
same as 6*3 in Eq.(6). The bubble grows beyond the critical size when the latent heat is conducted away from the 
surface into the surrounding medium which is governed by thermal dissipation and viscous damping. For our case, 
in the general framework of transition from a hadronic system to the QGP phase, we will use the expression for the 
dynamical prefactor from ref. [26] 

Here a is the surface tension of the bubble wall, Aw is the difference in the enthalpy densities of the QGP and the 
hadronic phases, A is thermal conductivity. Re is the critical bubble radius, and 77 and ^ are shear and bulk viscosities. 
(^ will be neglected as it is much smaller than ry. For A and 77, the following parametrizations are used [261 I27j . 

"^'[T^jynoJ ^\1 + Tyi000j\nj + 1 + 160/r2 ^ > 

/ai5\ /^y-^ / 0.02 \ { n^-^ 0.0225ri/2 

Here n/no is the ratio of the baryon density of the system to the normal nuclear baryon density, T is in MeV, ij 
is in MeV/fm^c, and A is in c/fm^. With this, the rate in Eq.(lO) is in fm~'*. For the range of temperatures of our 
interest (T > 160MeV), and for the low baryon density central rapidity region under consideration, it is the last n 
independent term for both 77 and k which dominates, and we will use these terms only for calculating 77 and k for our 
case. 

For the statistical prefactor, we use the following expression [35] 

The correlation length in the hadronic phase, £,hadi is expected to be of order of 1 fm and we will take it to be 0.7 
£tn [2S]). 

We will present estimates of the nucleation rates from Eq.(6) as well as Eq.(lO). One needs to determine the critical 
bubble profile and its 3-dimensional Euclidean action S^ (equivalently, AF in Eq.(lO)). We solve Eq.(7) using the 
fourth order Runge-Kutta method with appropriate boundary conditions (Eq.(8)), to get the profile of critical bubble 
[3T]. The critical bubble profiles (for the 3+1 dimensional case) are shown in Fig. (2a) for different temperatures. The 
bubble size decreases as temperature increases, since the energy difference between true vacuum and false vacuum 
increases (relative to the barrier height) as temperature increases. We choose a definite temperature T — 200 MeV for 
the nucleation of bubbles, which is suitably away from Tc to give acceptable bubble size and nucleation probabilities 
for the relevant time scale. Making T larger (up to T = 250 MeV when the barrier disappears) leads to similar bubble 
profile, and the nucleation probability of same order. 

Recall that we are calculating bubble profiles using Eq.(7) relevant for the 3+1 dimensional case, however, the field 
evolution is done using 2+1 dimensional equations as appropriate for the mid rapidity region of rapidly longitudinally 
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FIG. 2: 

(a) Critical bubble profiles for diflferent values of the temperature, (b) Solid curve shows the critical bubble for the 3+1 

dimensional case (which for finite temperature case becomes 3 Euclidean dimensional) for T = 200 MeV and the dotted curve 

shows the same for 2+1 dimensional case (i.e. 2 Euclidean dimensions for finite temperature case). 



expanding plasma. In Fig. (2b) solid curve shows the critical bubble for the 3+1 dimensional case (which for finite 
temperature case becomes 3 Euclidean dimensional) for T = 200 MeV and the dashed curve shows the same for 
2+1 dimensional case (i.e. 2 Euclidean dimensions for finite temperature case). It is clear that the 3 dimensional 
bubble is of supercritical size and will expand when evolved with 2+1 dimensional equations. This avoids the artificial 
construction of suitable supercritical bubbles which can expand and coalesce as was done in ref.|21J. (Recall that 
for the finite temperature case, a bubble of exact critical size will remain static when evolved by the field equations. 
In a phase transition, bubbles with somewhat larger size than the critical size expand, while those with smaller size 
contract.) 

For the bubble profile given by the sohd curve in Fig. (2b), the value of the action 83(1) is about 240 MeV. Using 
Eq.(6) for the nucleation rate, we find that the nucleation rate of QGP bubbles per unit time per unit volume is of 
the order of 0.025/?7i~''. The thermalization time for QGP phase is of the order of 1 fm at RHIC (say, for Au-Au 
collision at 200 GeV energy). Hence the time available for the nucleation of QGP bubbles is at most about 1 fm. We 
take the region of bubble nucleation to be of thick disk shape with the radius of the disk (in the transverse direction) 
of about 8 fm and the thickness of the disk (in the longitudinal direction) of about 1 fm. Total space-time volume 
available for bubble nucleation is then about 200 fm'' (in practice, less than this). For the case of Eq.(6), net number 
of bubbles is then equal to 5. 

For the case of the nucleation rate given by Eq.(lO), one needs an estimate of the critical bubble size Re as well 
as bubble surface energy a, (along with other quantities like A etc. as given by Eqs.(ll)-(14)). Determination of Re 
is somewhat ambiguous here as the relevant bubbles are thick wall bubbles as seen in Fig. (2). Here there is no clear 
demarcation between the core region and the surface region which could give an estimate of Re- Essentially, there is 
no core at all and the whole bubble is characterized by the overlap of bubble wall region. We can take, as an estimate 
for the bubble radius Re any value from 1 fm - 1.5 fm. It is important to note here that this estimate of Re is only 
for the calculation of nucleation rate F, and not for using the bubble profile for actual simulation. When bubbles are 
nucleated in the background of false vacuum with Z = 0, a reasonably larger size of the bubble is used so that cutting 
off the profile at that radius does not lead to computational errors and field evolution remains smooth. 

Once we have an estimate of Re, we can then estimate the surface tension a (which also is not unambiguous here) 



10 

as follows. With the realization that essentially there is no core region for the bubbles in Fig(2), we say that the 
entire energy of the bubble (i.e. the value of 53) conies from the surface energy. Then we write 

A-KRla = S3 (15) 

For S3 = 240 MeV, we get a = 8 MeV/fm^ if we take Re = 1.5 fm. With R^ = 1.0 fm, we get cr = 20 MeV/fm^. 

Number of bubbles expected can now be calculated for the case when nucleation rate is given by Eq.(lO). We find 
number of bubbles to be about lO^'* with Re taken as 1.5 fm. This is in accordance with the results discussed in ref. 
[26j . The bubble number increases by about a factor 5 if Re is taken to be about 1 fm. Thus, with the estimates 
based on Eq.(lO), bubble nucleation is a rare event for the time available for RHICE. 

As we have mentioned above, for us the bubble nucleation on one hand represents the possibility of actual dynamics 
of a first order transition, while on the other hand, it represents the generic properties of the domain structure arising 
from a C-D transition, which may very well be a cross-over, occurring in a finite time. With this view, and with 
various uncertainties in the determination of pre-exponential factors in the nucleation rate, we will consider a larger 
number of bubbles also and study domain wall and string production. First we will consider nucleation of 5 bubbles 
and then we will consider nucleation of 9 bubbles to get a better network of domain walls and strings. 

V. NUMERICAL TECHNIQUES 

In our simulation critical bubbles are nucleated at a time when the temperature T crosses the value T — 200 MeV 
during the initial stage between r = to r = tq = 1 fm (during which we have modeled the system temperature 
to increase linearly from to Tq). We take Tq = 400 MeV so bubble nucleation stage is taken to be at r = 0.5 fm 
when T reaches the value 200 MeV. Again, this is an approximation since in realistic case bubbles will nucleate over 
a span of time given by the (time dependent) nucleation rate, which could lead to a spectrum of sizes of expanding 
bubbles at a given time. However, due to very short time available to complete the nucleation of QGP bubbles in the 
background of confined phase, bubbles will have very little time to expand during the nucleation period of all bubbles 
(especially as initial bubble expansion velocity is zero) . Thus it is reasonable to assume that all the bubbles nucleate 
at the same time. 

After nucleation, bubbles are evolved by the time dependent equations of motion in the Minkowski space [21] as 
appropriate for Bjorken's longitudinal scaling model. 

.. /. d\ d\ -dV{l) 

with / = at T = 0. Here I = h + il2, and dot indicates derivative with respect to the proper time t. 

The bubble evolution was numerically implemented by a stabilized leapfrog algorithm of second order accuracy 
both in space and in time with the second order derivatives of li approximated by a diamond shaped grid. Here 
we follow the approach described in [3T] to simulate the first order transition. We need to nucleate several bubbles 
randomly choosing the corresponding Z(3) vacua for each bubble. This is done by randomly choosing the location of 
the center of each bubble with some specified probability per unit time per unit volume. Before nucleating a bubble, 
it is checked if the relevant region is in the false vacuum (i.e. it does not overlap with some other bubble already 
nucleated). In case there is an overlap, the nucleation of the new bubble is skipped. The orientation of I inside each 
bubble is taken to randomly vary between the three Z(3) vacua. 

For representing the situation of relativistic heavy ion collision experiments, the simulation of the phase transition is 
carried out by nucleating bubbles on a square lattice with physical size of 16 fm within a circular boundary (roughly 
the Gold nucleus size). We use fixed boundary condition, free boundary condition, as well as periodic boundary 
condition for the square lattice. To minimize the effects of boundary (reflections for fixed boundary, mirror reflections 
for periodic boundary conditions), we present results for free boundary conditions (for other cases, the qualitative 
aspects of our results remain unchanged). Even for free boundary conditions, spurious partial reflections occur, and 
to minimize these effects we use a thin strip (of 10 lattice points) near each boundary where extra dissipation is 
introduced. 

We use 2000 x 2000 lattice. For the physical size of 16 fm, we have Ax = 0.008 fm. To satisfy the Courant stability 
criteria, we use At = Ax/\^, as well as At = 0.9Aa;/-\/2, (which we use for the results presented in the paper). For 
Au-Au collision at 200 GeV, the thermalization is expected to happen within 1 fm time. As mentioned above, in this 
pre-equilibrium stage, we model the system as being in a quasi-cquilibrium stage with a temperature which increases 
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linearly with time (for siniplicity) . The temperature of the system is taken to reach upto 400 MeV in 1 fm time, 
starting from T = 0. After r = tq = 1 fm, the temperature decreases due to continued longitudinal expansion, i.e. 



T{t) = T{ro) (^) (17) 



r 

The stability of the simulation is checked by checking the variation of total energy of the system during the 
evolution. The energy fluctuation remains within few percent, with no net increase or decrease in the energy (for 
fixed and periodic boundary conditions, and without the dissipative I term in Eq.(16)) showing the stability of the 
simulation. 

The bubbles grow and eventually start coalescing, leading to a domain like structure. Domain walls are formed 
between regions corresponding to different Z(3) vacua, and strings form at junctions of Z(3) domain walls. Recall that 
the domain wall network is formed here in the transverse plane, appearing as curves. These are the cross-sections 
of the walls which are formed by elongation (stretching) of these curves in the longitudinal direction into sheets. At 
the intersection of these walls, strings form. In the transverse plane, these strings looks like vortices, which will be 
elongated into strings in the longitudinal direction. 

VI. EFFECTS OF QUARKS 

We now discuss the effects of quarks. As we mentioned above, we will follow the approach where the presence of 
quarks is interpreted as leading to explicit breaking of the Z(3) symmetry, lifting the degeneracy of different Z(3) 
vacua [5HH]- This has important effects in the context of our model. First of all, different vacua having different 
energies implies different nucleation rates for the QGP bubbles with different Z(3) vacua. Further, for non-degenerate 
vacua, even planar Z{3) interfaces do not remain static, and move away from the region with the unique true vacuum. 
Thus, while for the degenerate vacua case every closed domain wall collapses, for the non-degenerate case this is not 
true any more. A closed wall enclosing the true vacuum may expand if it is large enough so that the surface energy 
contribution does not dominate (this is essentially the same argument as given for the bubble expansion, see Eq.(5) 
and the discussion following it). 

To see the importance of these effects we need an estimate of the explicit symmetry breaking term arising from 
inclusion of quarks. For this we use the estimates given in O [25]. Even though the estimates in ref.[3nj are given 
in the high temperature limit, we will use these for temperatures relevant for our case, i.e. T ~ 200 — 400 MeV, to 
get some idea of the effects of the explicit symmetry breaking. The difference in the potential energy between the 
true vacuum with I = 1 and the other two vacua {I = e*^^/^, and / — e'^'^/'', which are degenerate with each other) is 
estimated in ref. [53] to be, 

AV^^7:'T^^iN'-2) (18) 

where Ni is the number of massless quarks. If we take Ni ^ 2 then AV" ~ 3T'*. At the bubble nucleation temperature 
(which we have taken to be about T = 200 MeV, the difference between the false vacuum and the true vacuum is 
about 150 MeV/fm'^ while AT^ at T = 200 MeV is about four times larger, equal to 600 MeV/fm'^. As T approaches 
Tc, this difference will become larger as the metastable vacuum and the stable vacuum become degenerate at Tc, while 
Ay remains non-zero. For T near 250 MeV (where the barrier between the metastable vacuum and the stable vacuum 
disappears), AV becomes almost comparable to the difference between the potential energy of the false vacuum (the 
confining vacuum) and the true vacuum (deconfined vacuum). 

It does not seem reasonable that at temperatures of order 200 MeV a QGP phase (with quarks) has higher free 
energy than the hadronic phase. This situation can be avoided if the estimates of Eq.(18) are lowered by about a 
factor of 5 so that these phases have lower free energy than the confining phase. (A more desirable situation will be 
when AV approaches zero as the confining vacuum and the deconfining (true) vacuum become degenerate at Tc.) It is 
in the spirit of the expectation that explicit breaking of Z(3) is small near Tc for finite pion mass [H]- Even with such 
lower estimates, the effects of quarks may give different nucleation probabilities for different Z(3) vacua. However, 
in this paper we will ignore this possibility. This may not be very unreasonable as for thick wall bubbles thermal 
fluctuations may be dominant in determining the small number of bubbles nucleated during the short span of time 
available. For very few bubbles nucleated, there may be a good fraction of events where different Z(3) vacua may 
occur in good fraction. Also note that the pre-exponential factor for the bubble nucleation rate of Eq.(6), as given 
in Eq.(9) increases with the value of S3{1). Thus, for the range of values of 83(1) for which the exponential factor in 
Eq.(6) is of order 1, which is likely in our case, the nucleation rate may not decrease with larger values of 83(1), i.e. 
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for the Z(3) vacua with higher potential energies than the true vacuum. (Of course for very large values of S3{1) the 
exponential term will suppress the nucleation rate.) Thus our assumption of neglecting quark effects for the bubble 
nucleation rate may not be unreasonable. 

We now consider the effect of non-zero AV, as in Eq.(18), on the evolution of closed domain walls. The temperature 
range relevant for our case is T = 200 — 400 McV. In an earlier work we had numerically estimated the surface tension 
of Z(3) walls to be about 0.34 and 7.0 GeV/fm^ for T = 200 and 400 MeV respectively. The effects of quarks wih be 
significant if a closed spherical wall (with true vacuum inside) starts expanding instead of collapsing. Again, using the 
bubble free energy Eq.(5), with 77 = AV and a as the surface energy of the interface, we see that the critical radius 
R* of the spherical wall is 

R* = — ^ 2a/3T^ (19) 

V 

For T = 200 MeV and 400 MeV we get i?* ~ 1 and 1.5 fm respectively. Though these values are not large, 
these are not too small either when considering the fact that relevant sizes and times for RHICE are of order few 
fm anyway. The values of R* we estimated here are very crude as for these sizes wall thickness is comparable to R* 
hence application of Eq.(5), separating volume and surface energy contributions, is not appropriate. Further, as we 
discussed above, the estimate in Eq.(18) which is applicable for high temperature limit, seems an over estimate by 
about an order of magnitude at these temperatures. Thus uncertainties of factors of order 1 may not be unreasonable 
to expect. In that case the dynamics of closed domain walls of even several fm diameter will not be affected by the 
effects of quarks via Eq.(18). 

We will see in the next section of simulation results that the domain walls and strings typically have large velocities 
(e.g. about 0.5 - 0.8) at the time of formation. These result from momentum of colliding bubble walls and from 
curvature in the shape of these walls (as well as asymmetries in the profiles of strings) at the time of formation. With 
such large velocities present, the effects of pressure differences between different Z(3) vacua due to quarks may become 
subdominant in studying the evolution of these structures for the short time duration available for RHICE. 

With this, we will assume that for small closed walls, of order few fm diameter, as is expected in RHICE, the quark 
effects in the evolution of wall network may be neglected. We plan to remove these assumptions and include the 
effects of non-zero AV^ due to quarks on bubble nucleation and wall evolution in a future work. 

VII. RESULTS OF THE SIMULATION 

As we mentioned above, the number of bubbles expected to form in RHICE is small. We first present and discuss 
the case of 5 bubbles which is more realistic from the point of view of nucleation estimates given by Eq.(6) (though a 
gross over-estimate for Eq.(lO)). Note that a domain wall will form even if only two bubbles nucleate (with different 
Z(3) vacua). However, to see QGP string formation, we need nucleation of at least three bubbles. Next we will discuss 
the case of 9 bubbles which is a much more optimistic estimate of the nucleation rate (even for Eq.(6)). Alternatively, 
this case can be taken as better representation of the case when the transition is a cross over and bubbles only 
represent a means for developing a domain structure expected after the cross-over is completed. (In this case only 
relevant energy density fluctuations, as discussed below, will be those arising from Z(3) walls and strings, and not the 
ones resulting from bubble wall coalescence.) 

A. formation and motion of extended walls 

Fig. 3-7 show the results of simulation when five bubbles are nucleated with random choices of different Z(3) vacua 
inside each bubble. Fig. 3 shows a time sequence of surface plots of the order parameter l{x) in the two dimensional 
lattice. Fig. 3a shows the initial profiles of the bubbles of the QGP phase embedded in the confining vacuum with 
Z = at r = 0.5 fm with the temperature T = 200 MeV. (Recall that for initial 1 fm time the temperature is taken 
to linearly increase from zero to Tq = 400 MeV.) The radial profile of each bubble is truncated with appropriate 
care of smoothness on the lattice for proper time evolution. Fig. 3b shows the profile of each bubble at r = 1.5 fm 
showing the expansion of the bubbles. Near the outer region of a bubble the field grows more quickly towards the 
true vacuum. If bubbles expand for long time then bubble walls become ultra-relativistic and undergo large Lorentz 
contraction. This causes problem in simulation (see, e.g. |21j). In our case this situation arises at outer boundaries 
(for the inner regions bubbles collide quickly) . For outer regions also it does not cause serious problem because of the 
use of dissipative boundary strip (as explained above in Sect.V). 
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FIG. 3: 

(a) and (b) show plots of profiles oi I at t — 0.5 fm and 1.5 fm respectively, (c)-(i) show plots oi lo — I at t = 1.5, 2.5, 4.0, 6.0, 

9.0, 11.0, and 13.7 fm. T drops to below Tc around at r = 10.5 fm and T = 167 MeV at r = 13.7 fm. Formation of domain 

walls and string and antistring (at junctions of three walls) can be seen in the plots in (e) - (h). 



Figs. 3c-3i show plots of ^o — ^ clearly showing formation of domain walls and strings (junctions of three walls). 
Here Iq is a reference vacuum expectation value of I calculated at the maximum temperature T = Tq = 400 MeV. 
Formation of domain walls, extending through the entire QGP region, is directly visible from Fig.Se (at t = 4 fm) 
onwards. The temperature drops to below Tc ~ 182 MeV at t ~ 10.5 fm. The last plot in Fig.3i is at r = 13.7 fm 
when the temperature T = 167 MeV, clearly showing that the domain walls have decayed away in the confined phase 
and the field is fiuctuating about I ~ 0. 

In Fig.3f-3h we see two junctions of three domain walls where the QGP strings form. This is seen more clearly in 
Fig. 4 where the phase 6 oi I is plotted (with the convention that 9 is the angle of the arrow from the positive X axis). 
The domain walls are identified as the boundaries where two different values of 9 meet, and strings correspond to 
the non-trivial winding of 9 at the junctions of three walls. From Fig. 4b, c we clearly see that at one of the junctions 
we have a string (at X ~ 5 fm, K ~ 8 fm) with positive winding and we have an anti-string at X ~ 9 fm, 1" ~ 8 
fm with negative winding. Note the rapid motion of the walls forming the antistring towards positive Y axis from 
Fig. 4c (at r = 4.6 fm) to Fig.4e (at r = 12.2 fm). The average speed of the antistring (and wall associated with that) 
can be directly estimated from these figures to be about 0.5 (in natural units with c = 1). This result is important 
in view of the discussion in the preceding section showing that effects of pressure differences between different Z(3) 
vacua, arising from quarks, may be dominated by such random velocities present for the walls and strings at the time 
of formation. The motion of the walls here is a direct result of the straightening of the L shaped wall structure due 
to its surface tension. For the same reason the wall in the left part of Fig. 4b also straightens from the initial wedge 
shape. Fig.4f shows almost random variations oi 9 at t — 13.7 fm when the temperature is 167 MeV, well below the 
critical temperature. Though it is interesting to note that a large region of roughly uniform values of 9 still survives 
at this stage. 

Fig. 5 shows surface plots of the local energy density e. e is plotted in units of GcV/fm'^. Although the simulation 
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FIG. 4: 

Plots of the phase 6 of the order parameter /. (a) shows the initial distribution of 6 in the bubbles ai t — 0.5 fm. (b) - (f) 

show plots of S at T = 2.0, 4.6, 11.0, 12.2, and 13.7 fm respectively. Location of domain walls and the string (with positive 

winding) and antistring (with negative winding) are clearly seen in the plots in (b)-(e). The motion of the antistring and 

associated walls can be directly seen from these plots and an estimate of the velocity can be obtained. 



is 2+1 dimensional representing the transverse plane of the QGP system, we calculate energy in 3+1 dimensions by 
taking a thickness of 1 fm in the central rapidity region. Fig. 5a shows plot at r = 3 fm when bubbles have coalesced. 
In Fig. 5b at t = 3.6 fm we see that bubble walls have almost decayed (in ripples of I waves) between the bubbles 
with same 9 (i.e. same Z(3) vacua) as can be checked from 9 plots in Fig. 4. Energy density remains well localized 
in the regions where domain walls exist. Also, one can see the small peaks in the energy density where strings and 
antistrings exist. Large peaks arise from oscillations of I when bubble walls coalesce, as discussed in [21,. Large values 
of e near the boundary of the lattice are due to relativistically expanding bubble walls. Motion of walls and generation 
of increased fluctuations in energy density are seen in Fig.5c-5e. Fig.5f at r = 13.2 fm (with T — 169 MeV) shows 
that walls have decayed. However, some extended regions of high energy density can be seen at this stage also. 

Fig. 6 shows contour plots of energy density e. Fig. 6a shows coalescence of bubbles at r = 3 fm. Fig. 6b at t — 
3.6 fm clearly shows the difference in the wall coalescence depending on the vacua in the colliding bubbles. Where 
domain walls exist, we see extended regions of high energy density contours whereas where the two vacua in colliding 
bubbles are same, there are essentially no high energy density contours. Motion of domain walls (and strings at wall 
junctions) is clearly seen in these contour plots in Fig.6b-6d. Fig.6e is at t = 10.2 fm when T drops to Tc- Wall 
structures are still present. Fig.6f is at t = 13.2 fm [T = 169 MeV) when walls have decayed away, though some 
extended structures in contours still survive. 

We have also calculated the variance of energy density Ae at each time stage to study how energy fluctuations 
change during the evolution. In Fig. 7 we show the plot of Ae/e as a function of proper time. Here e is the average 
value of energy density at that time stage. The energy density e decreases due to longitudinal expansion, hence we 
plot this ratio to get an idea of relative importance of energy density fluctuations. Fig. 7 shows initial rapid drop in 
Ae/e due to large increase in e during the heating stage upto t = 1 fm, followed by a rise due to increased energy 
density fluctuations during the stage when bubbles coalesce and bubble walls decay, as expected. Interesting thing 
to note is a slight peak in the plot near r = 10.5 fm when T drops below Tc- This should correspond to the decay of 
domain walls and may provide a signal for the formation and subsequent decay of such objects in RHICE. 
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FIG. 5: 

Surface plots of the local energy density e in GeV/fm^. (a) - (f) show plots at r = 3.0, 3.6, 5.0, 6.0, 8.0, and 13.2 fm 

respectively. Extended domain walls can be seen from these plots of e in (b) - (e). Small peaks in e exist at the locations of 

string and antistring (larger peaks arise from oscillations of field where bubbles coalesce). Plot in (f) is at the stage when T = 

169 MeV and domain walls have decayed away. 



B. Formation and collapse of a closed wall 



Fig. 8-12 show the results of simulation where nine bubbles are nucleated. Fig. 8 shows a time sequence of surface 
plots of l{x) (similar to Fig. 3). Fig. 8a shows the initial profiles of the QGP bubbles at r = 0.5 fm with the nucleation 
temperature of T = 200 MeV. Fig. 8b shows the profile of I for the bubbles at r = 1.5 fm showing the expansion of 
the bubbles. Fig.8c-8i show plots oi Iq — I at different stages. Noteworthy here is the formation of a closed domain 
wall near the central region which is clearly first seen in Fig.8e at r = 5 fm. The collapse of this closed domain wall is 
seen in the subsequent plots with the closed wall completely collapsing away in Fig.8h at r = 9.6 fm. Only surviving 
structure is an extended domain wall along the X axis. Fig.8i is at r = 13.2 fm when T = 169 MeV. The domain 
walls have decayed away and / fluctuates about the value zero as appropriate for the confined phase. 

Fig. 9 shows plots of the phase of / at different stages. Initial phase distribution in different bubbles is shown in 
Fig. 9a at r = 0.5 fm. Fig. 9b shows the formation of closed, elliptical shaped domain wall at r = 2.6 fm. Strings 
and antistrings can also be identified by checking the windings of 0. The closed domain wall collapses, and in the 
process becomes more circular, as shown in the plots in Fig.9b-9h. Fig.9h shows the plot at r = 9.6 fm when the 
closed domain wall completely collapses away, leaving only an extended domain wall running along X axis between 
y ~ 4 — 8 fm. The final Fig.9i si t — 13.8 fm is when the temperature T — 167 MeV showing random fluctuations 
of 9 when domain walls have decayed away in the conflning phase. 

Fig. 10 shows the surface plot of energy density e at different stages. Extended thin regions of large values of e are 
clearly seen in the plots corresponding to domain walls. Collapse of the closed domain wall is also clearly seen in 
Fig.lOc-lOg. Important thing to note here is the surviving peak in the energy density plot at the location of domain 
wall collapse. This peak survives even at the stage shown in Fig.lOi at r = 13.2 fm when T = 169 MeV, well below 
the transition temperature. Such hot spots may be the clearest signals of formation and collapse of Z(3) walls. 

Contour plots of e are shown in Fig. 11. Though closed domain wall can be seen already in Fig. lib (at r — 3.5 
fm), the domain wall is still attached to outward expanding bubble walls near X = i fm, y = 12 fm which affects 
the evolution/motion of that portion of the domain wall. Formation of distinct closed wall structure is flrst visible in 
Fig. lie at T =4.6 fm. Subsequent plots clearly show how the domain wall becomes circular and finally collapses away 
by Fig. Ilk at r = 9.6 fm. Note the survival of the hot spot even at the stage shown in Fig. 11/ at r = 13.2 fm when 
T = 169 MeV. 
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FIG. 6: 
Contour plots of the local energy density e at different stages. Plots in (a) - (f) correspond to r = 3.0, 3.6, 5.0, 7.6, 10.2, and 
13.2 fm respectively. Structure of domains walls formed near the coalescence region of bubbles with different Q is clear in (b), 
whereas the bubble walls at lower half of Y region, and near X = 10 fm, are seen to simply decay away due to same vacuum 
in the colliding bubbles. Motion of the antistring and associated domain walls is clear from plots in (b) - (e). The last plot in 

(f) is at r = 13.2 fm when T = 169 MeV. 
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FIG. 7: 

Plot of the ratio of variance of energy density Ae and the average energy density e as a function of proper time. Energy 

fluctuations increase during the initial stages when bubbles coalesce and bubble walls decay. After that there is a slow 

decrease in energy fluctuation until the stage when the temperature drops below Tc and r ~ 10.5 fm. Energy fluctuations 

increase after this stage. Note small peak near the transition stage. 
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FIG. 8: 
(a) and (b) show plots of profiles of Z at r = 0.5 fm and 1.5 fm respectively for the case when 9 bubbles are nucleated, (c)-(i) 

show plots of Zo — ^ at r = 2.0, 3.0, 5.0, 7.0, 8.6, 9.6, and 13.2 fm respectively. Formation of a closed domain wall is first 

clearly seen in the plot in (e). This closed domain wall collapses as seen in plots in (e) through (h). Only surviving domain 

wall is an extended wall along X axis in (h). Plot in (i) is when the temperature T = 169 MeV. 



One can make a rough estimate of the velocity of the closed wall during its collapse from these plots. In Fig. lie, 
aX T — 4.6 fm, the X extent of the closed wall is about 8 fm and the Y extent is about 5 fm. The wall collapses 
away by the stage in Fig. Ilk at r = 9.6 fm. This gives rough velocity of collapse in X direction to be about 0.8 while 
the velocity in Y direction is about 0.5. Note that here, as well as in Fig. 4 for the five bubble case, the estimate of 
the wall velocity is not affected by the extra dissipation which is introduced only in a very thin strip (consisting of 
ten lattice points) near the lattice boundary. Formation and collapse of such closed domain walls is important as the 
resulting hot spot can lead to important experimental signatures. Further, such closed domain wall structures are 
crucial in the studies of Pt enhancement, especially for heavy flavor hadrons as discussed in [TTJ [T^] 

Fig. 12 shows the evolution of the ratio of the variance of energy density and the average energy density. As for the 
five bubble case, initial drop and rise are due to heating stage upto t = 1 fm and subsequent bubble coalescence and 
decay of bubble walls. In this case the ratio remain roughly constant upto t ~ 10.5 fm which is the transition stage 
to the confining phase. This is the stage when the surviving extended domain wall starts decaying. This is also the 
stage soon after the closed domain wall collapses away. The prominent peak at this stage should be a combined result 
of both of these effects. The large increase in the variance of energy density at this stage should be detectable from 
the analysis of particle distributions and should be a clear signal of hot spots resulting from collapse of closed walls 
and the decay of any surviving domain walls. 
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FIG. 9: 

Plots of the phase Q of /. (a) shows initial distribution of Q in the bubbles at r = 0.5 fm. (b) - (i) show plots of 61 at r = 2.6, 

5.0, 6.0, 7.0, 8.0, 8.6, 9.6, and 13.8 fm respectively, (b) shows formation of elliptical shaped closed domain wall which 

subsequently becomes more circular as it collapses away by r = 9.6 fm as shown by the plot in (h). The plot in (i) is at 

T = 167 MeV showing random fluctuations of Q. 

VIII. POSSIBLE EXPERIMENTAL SIGNATURES OF Z(3) WALLS AND STRINGS 



The Z(3) wall network and associated strings form during the early confincmcnt-deconfinenient phase transition. 
They undergo evolution in an expanding plasma with decreasing temperature, and eventually melt away when the 
temperature drops below the deconfinement-confinement phase transition temperature. They may leave their signa- 
tures in the distribution of final particles due to large concentration of energy density in extended regions as well as 
due to non-trivial scatterings of quarks and antiquarks with these objects. 

First, we focus on the extended regions of high energy density resulting from the domain walls and strings. This 
is clearly seen in our simulations and some extended structures/hot spots also survive after the temperature drops 
below the transition temperature T^. Note that even the hot spot resulting from the collapse of closed domain wall in 
Fig. 9, 10 will be stretched in the longitudinal direction into an extended linear structure (resulting from the collapse of 
a cylindrical wall). We know that at RHIC energies, the final freezeout temperature is not too far below the transition 
temperature Tj,. This means that the energy density concentrated in any extended (sheet like for domain walls and 
line like for strings/hot spots) regions may not be able to defuse away effectively. Assuming local energy density to 
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FIG. 10: 

Surface plots of the local energy density e in GeV/fm'^. (a) - (i) show plots at r = 2.6, 3.5, 4.6, 5.6, 7.0, 8.6, 9.6, 11.2, and 

13.2 fm respectively. Formation and subsequent collapse of closed domain wall is clearly seen in plots in (c) through (g). Note 

that the strong peak in e resulting from domain wall collapse (the hot spot) survives in (i) when T — 169 MeV. 



directly result in multiplicity of particles coming from that region, an analysis of particle distribution in Pt and in 
rapidity should be able to reflect any such extended regions. In this context, it will be interesting to investigate if the 
ridge phenomenon seen at RHIC [3D] could be a manifestation of an underlying Z(3) domain wall/string structure. 
Correlation of particle production over large range of rapidity will naturally result from longitudinally extended 
regions of high energy density (hot spots in the transverse plane). Combined with flow effects it may lead to ridge like 
structures [301 131) . If extended domain wall structure survives in the transverse plane also, this will then extend to 
sheet like regions in the longitudinal direction. Decay of such a region of high energy density may directly lead to a 
ridge like structure, without requiring flow effects. 

We expect non-trivial signatures resulting from the consideration of interactions of quarks and antiquarks with 
domain walls. It was shown in an earlier work [TT] using generic arguments that quarks and antiquarks should 
have non-zero reflection coefficients when traversing across these domain walls. A collapsing domain wall will then 
concentrate any excess baryon number enclosed, leading to formation of baryon rich regions. This is just like Witten's 
scenario for the early universe ^^ (which was applied for the case of RHICE in ref. [33j). However, for these works 
it was crucial that the quark-hadron transition be of (strong) first order. As we have emphasized above, in our case 
formation of Z(3) walls and strings will be a generic feature of any C-D phase transition. Even though we have 
implemented it in the context of a first order transition via bubble nucleation, these objects will form even if the 
transition is a cross-over. Thus, concentration of baryons in small regions should be expected to occur in RHICE 
which should manifest in baryon concentration in small regions of rapidity and Pt ■ 

Another important aspect of quark/antiquark reflection is that inside a collapsing wall, each reflection increases 
the momentum of the enclosed particle. When closed domain walls collapse then enclosed quarks/antiquarks may 
undergo multiple reflections before finally getting out. This leads to a specific pattern of Pt enhancement of quarks 
with heavy flavors showing more prominent effects |12) . The modification of Pt spectrum of resulting hadrons can be 
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FIG. 11: 
Contour plots of the local energy density e at different stages. Plots in (a) - (1) correspond to r = 2.6, 3.5, 4.6, 6.0, 6.6, 7.0, 
7.6, 8.0, 8.6, 9.0, 9.6, and 13.2 fm respectively. Formation of distinct closed wall structure is first visible in (c) at r =4.6 fm. 
Subsequent plots show the collapse of this domain wall as it becomes more circular. The wall finally collapses away in (k) at 
T = 9.6 fm. Note that concentration of energy density at the location of domain wall collapse (the hot spot) survives even at 

the stage shown in (1) at r = 13.2 fm when T = 169 MeV. 
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FIG. 12: 

Plot of the ratio of variance of energy density Ae and the average energy density £ as a function of proper time. Energy 
fluctuations increase during the initial stages when bubbles coalesce and bubble walls decay. After that Ae/e remains roughly 
constant until the stage when the temperature drops below Tc at r ~ 10.5 fm. This is also the stage just after the collapse of 

the closed domain wall. Energy fluctuations sharply increase around this stage. Note the prominent peak at this stage. 

calculated, and the enhancement of heavy flavor hadrons at high Pt can be analyzed for the signal for the formation 
of Z(3) domain walls in these experiments [12!. In our simulations extended domain walls also form which show 
bulk motion with velocities of order 0.5. Quarks/antiquarks reflected from such moving extended walls will lead to 
anisotropic momentum distribution of emitted particles which may also provide signature of such walls. For collapsing 
closed domain walls, spherical domain walls were used for estimates in ref. |llj and in ref. |12j . Our simulation in the 
present work provides a more realistic distribution of shapes and sizes for the resulting domain wall network. We have 
estimated the velocity of moving domain walls to range from 0.5 to 0.8 for the situations studied. These velocities 
are large enough to have important effect on the momentum of quarks/antiquarks undergoing reflection from these 
walls. One needs to combine the analysis of [TT1[T2] with the present simulation to get a concrete signature for baryon 
concentration and heavy flavor hadron Pt spectrum modification. We plan to carry this out in a future work. We 
also plan to study effects of spontaneous violation of CP due to formation of these Z(3) walls in RHICE. 

Our results show interesting pattern of the evolution of the fluctuations in the energy density. As seen in Fig. 7 
and Fig. 12, energy density fluctuations show rapid changes during stages of bubble wall coalescence and during 
collapse/decay of domain walls. Even string-antistring annihilations should be contributing to these fluctuations. 
Fluctuations near the transition stage may leave direct imprints on particle distributions. It is intriguing to think 
whether dileptons or direct photons may be sensitive to these fluctuations, which could then give a time history of 
evolution of such energy density fluctuations during the early stages as well. Even the presence of domain walls and 
strings during early stages may affect quark-antiquark distributions in those regions which may leave imprints on 
dileptons/direct photons. An important point to note is that in our model, we expect energy density fluctuations 
in event averages (representing high energy density regions of domain walls/strings as discussed above), as well as 
event-by-event fluctuations. These will result due to fluctuation in the number/geometry of domain walls/strings 
from one event to the other resulting from different distribution of (randomly occurring) Z(3) vacua in the QGP 
bubbles. Even the number of QGP bubbles, governed by the nucleation probability, will vary from one event to the 
other contributing to these event-by-event fluctuations. 



IX. CONCLUSIONS 



We have carried out numerical simulation of formation of Z{3) interfaces and associated strings at the initial 
confinement-deconfinement phase transition during the pre-equilibrium stage in relativistic heavy-ion collision exper- 
iments. A simple model of quasi-equilibrium system was assumed for this stage with an effective temperature which 
first rises (with rapid particle production) to a maximum temperature Tq > Tc, and then decreases due to continued 
plasma expansion. 

Using the effective potential for the Polyakov loop expectation value l{x) from ref. [7|[H] we study the dynamics of 
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the (C-D) phase transition in the temperature/time range when the first order transition of this model proceeds via 
bubble nucleation. As we have emphasized above, though our study is in the context of a first order transition, its 
results are expected to be valid even when the transition is a cross-over. (Though for non-zero chemical potential the 
transition may indeed be of first order) . The generic nature of our results arises due to the fact that the formation of 
Z(3) domain walls and associated strings happens due to the a general domain structure resulting after any transition 
(occurring in a finite time). This is the essential physics of the Kibble mechanism underlying the formation of 
topological defects in symmetry breaking transitions. 

The Z{3) wall network and associated strings formed during this early C-D transition are evolved using field 
equations in a plasma which is longitudinally expanding, with decreasing temperature. We have neglected here the 
transverse expansion which is a good approximation for the early stages near the formation stage of these objects, 
but may not be a good approximation for the later parts of simulations when temperature drops below T^ and Z(3) 
domain walls and strings melt away. We have studied size/shape of resulting closed domain wall as well as extended 
domain walls and have estimated the velocities of walls to range from 0.5 to 0.8. We also calculate the energy density 
fluctuations expected due to formation of these objects. Various experimental signals which can indicate the formation 
of these topologically non-trivial objects in RHICE have been discussed. For example, existence of these objects will 
result in specific patterns of energy density fluctuations which may leave direct imprints on particle distributions. 
In our model, we expect energy density fluctuations in event averages (representing high energy density regions of 
domain walls/strings), as well as event-by-event fluctuations as the number/geometry of domain walls/strings and 
even the number of QGP bubbles, varies from one event to the other. Extended regions of large energy densities 
arising from Z(3) walls and associated strings may be manifested in space-time reconstruction of hadron density 
(using hydrodynamic model). The correlation of particle production over large range of rapidity will be expected 
from such extended regions. This, combined with the flow effects (for string like regions), or possibly directly (for 
sheet like extended region) may provide an explanation for the ridge phenomena observed at RHIC |30j . Also, from 
the reflection of quarks and antiquarks from collapsing domain walls, baryon number enhancement in localized regions 
(due to concentration of net baryon number) as well as enhancement of heavy flavor hadrons at high Pt is expected. 

We emphasize again that the presence of Z{3) walls and string may not only provide qualitatively new signatures 
for the QGP phase in these experiments, it may provide the first (and may be the only possible) laboratory study of 
such topological objects in a relativistic quantum field theory system. 
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